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Abstract. A two-dimensional generalized hydrodynamics (GH) model is developed to study 
the full spectrum of both longitudinal and transverse dust acoustic waves (DAW) in strongly 
coupled complex (dusty) plasmas, with memory-function-formalism being implemented to 
enforce high-frequency sum rules. Results are compared with earlier theories (such as quasi- 
localized charge approximation and its extended version) and with a self-consistent Brownian 
dynamics simulation. It is found that the GH approach provides good account, not only for 
dispersion relations, but also for damping rates of the DAW modes in a wide range of coupling 
strengths, an issue hitherto not fully addressed for dusty plasmas. 
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1. Introduction 

Dusty, or complex plasmas have grown into a mature research field with a surprisingly 
broad range of interdisciplinary facets E |2l O. There has been significant interest in the 
study of dynamics of collective excitations of dusty plasmas, both in experiments S [51 [H 
111 [H m [TOl and through theory/numerical simulation of strongly-coupled Yukawa systems 
|[IIl[I2[Il[Il[l5l[ia[ni[Il[Bl2ai2Dl22l|2l. Of particular interest in these studies 
is theoretical modeling of longitudinal and transverse dust acoustic wave (DAW) modes 
ifm . and consequently several assessments of how the existing analytical theories compare 
with the results from computer experiments have appeared recently lfT9l l20l l22l l23l l24l| . 
Specifically, by comparing the Brownian dynamics (BD) simulation results for the DAW 
spectra with various theoretical models, it was shown that the so-called Quasi-localized charge 
approximation (QLCA) ifTSl l20ll provides an overall good account of the DAW dispersion 
relations, i.e., the peak locations in those spectra, for both three-dimensional (3D) and two- 
dimensional (2D) dusty plasma states characterized by the coupling strength V in the range 
10 < r < 1000 [|22ll24l . However, since the QLCA provides no information about the profile 
of such spectra and, consequently, gives no account of damping processes for the collective 
modes [fT5l . it is desirable to explore alternate theoretical approaches to strongly coupled 
systems. In that context, generalized hydrodynamics (GH) offers possibility to describe 
full wave spectra by introducing viscoelastic effects in the dynamics of dusty plasmas in 
a phenomenological manner, as was demonstrated in studying DAWs in 3D dusty plasmas 

ttiaiiiini. 

In contrast to ordinary hydrodynamics (OH), which is only valid for weakly coupled 
fluids with r ^ 1 in the long wavelength limit, GH covers wider ranges of coupling 
strengths and wavelengths, where structural effects begin to show up. Particularly suitable and 
physically intuitive framework for implementing the GH model is provided by descriptions of 
the Lennard- Jones fluids, pioneered by Ailawadi et al. [|231 . Boon and Yip [|26l . and Hansen 
and McDonald fi27il . This framework may be applied to a dusty plasma by asserting that dust 
particles form a quasi-neutral fluid with the inter-particle interactions described via a Debye- 
Hiickel, or Yukawa potential, which arises due to screening of the charge accumulated on dust 
particles by the background electron and ion fluids. Successful applications of the GH model 
to dusty plasma were conducted by Kaw and Sen [fT3l and Murillo lfT4l to study the wave 
dispersion and the shear wave cutoff in 3D dust liquids, respectively. 

However, there is still great demand for a GH model for 2D strongly coupled dusty 
plasmas (SCDPs), especially because such systems have become particularly favored in recent 
laboratory experiments [|3l|28|. Therefore, our principal motivation in this work is to develop 
and implement a GH model for studying the collective dynamics in 2D SCDPs. In addition, 
it is necessary to generalize the GH model to include the effects of collisions of dust particles 
with the neutral gas in the background plasma in a manner similar to that used for colloidal 
suspensions in describing the collisions of macroions with the solvent molecules [291 . While 
such collisions give rise to the classical Brownian motion of macroparticles in both these 
systems, it should be stressed that the friction forces on dust particles due to the neutrals in 
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dusty plasmas are generally much weaker than the friction forces on the colloidal macroions. 
Accordingly, effects of the neutral gas on the DAW dispersion relations are found to be 
negligible for the range of parameters of interest in experiments with dusty plasmas, thus 
validating the use of molecular dynamics (MD) simulations for such systems [24J. However, 
it is still unclear as to what extent can damping due to the neutral gas compete with the viscous 
damping of longitudinal DAWs at finite frequencies and finite wavelengths in dusty plasmas 
[[30l. While this difficult issue naturally arises in this work, our primary motivation in using 
BD simulation lies in the fact that treating the dust particles as Brownian particles provides 
a natural and convenient way to eliminate the need for thermostation that arises in the MD 
simulation of dusty plasmas [[2211311 . 

In this paper we perform a BD simulation of strongly coupled 2D Yukawa liquids, 
and use the results to evaluate the equilibrium radial distribution function (RDF) and the 
static structure factor, as well as the (power) spectral densities for both the longitudinal and 
transverse current densities in such systems. The former two quantities enable us to use 
the GH approach within the memory function formalism to evaluate both the dispersion 
relation and the damping rate of the DAW modes in 2D dusty plasmas by enforcing the 
low-order, high-frequency sum rules upon the theoretical spectral densities. It is found that 
the GH approach provides dispersion relations that compare well with those resulting from 
the simulation spectra over broad ranges of wavelengths and coupling strengths. Additional 
comparison with the dispersion relations from the QLCA model shows that the GH approach 
provides a good account of the direct thermal effect at lower coupling strengths and shorter 
wavelengths, which is absent in the QLCA approach, but is seen in the simulation data. On 
the other hand, a simple extension of the QLCA dispersion relation, which was recently 
proposed in Ref. [[22|[. is found to be in a much better agreement with both the GH results 
and simulation data. Most importantly, the GH results are shown to yield a good fit to the 
spectral density profiles from our BD simulations, thus providing semi-analytical modeling 
of the wave-number dependent damping rates of the DAW modes in strongly coupled dusty 
plasmas in a broad range of coupling strengths, an issue that has not been fully addressed 
so far [[22|. Finally, we also find that the effect of collisions with the neutrals on the DAW 
damping rates is well accounted for by introducing a local in time friction force into the GH 
equations. 

The manuscript is organized as follows. The details of the BD simulation are given in 
Section 2. Theoretical development of the generalized hydrodynamics is discussed in Section 
3. Results based on two simple models for the memory function are presented in Section 4, 
while we conclude in Section 5. 

2. Simulation 

We perform BD simulation of a 2D system consisting of = 4000 dust particles, each 
carrying a constant charge g^, which are initially placed at random positions within a square 
simulation cell in the r = {x,y) plane with periodic boundary conditions. Evolution of the 
system is governed by equations that may be regarded as a stochastic generalization of the 
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usual MD method, where the effect of collisions with the neutral particles in the background 
plasma is modeled by a Langevin generalization of the usual Newton's equations. Therefore, 
equations for the velocity and the position vectors of the i-th dust particle are given by 

= - 7nv.(t) + ^F(r,(t)-r,(t)) + A,(t), 

where is the mass of each dust particle and F(r) = —{r/r)dU{r)/dr is the pair-wise 
force between two dust particles a distance r = a/^;^ + apart, interacting via the Yukawa 
potential, U{r) = (q'^/r) exp (— r/ A/)), with A/j being the Debye screening length of electrons 
and ions in the background plasma. Collisions of the i-th dust particle with the neutrals in the 
plasma give rise to a systematic drag force, which we model by the Epstein drag coefficient 
7„, and to a cumulative random force, which we model by a delta-correlated Gaussian white 
noise. When the system is in thermal equilibrium, the friction coefficient 7„ and the stochastic 
acceleration of the i-th particle Ai(t) are related to the ambient temperature Td via the 
fluctuation-dissipation theorem. Equations ^ are solved for our system by the Gear-like 
predictor-corrector algorithm for BD simulation [l3Tl, which was previously used in modeling 
both the equilibrium structure ||32l and collective modes in 2D dusty plasmas [[22l|. 

Equations ([T]) can be expressed in terms of just three parameters ll22l : the coupling 
strength F = g^/ (a/csTd), the reduced neutral drag coefficient 7 = jn/^pd, and the screening 
parameter k = a/ Ad, where a = 1/ y/Wp is the average inter-particle separation with p being 
the average surface density of dust particles, and cOpd = \/'2qj/ {nida^) is the characteristic 
dust-plasma frequency. In this work, we perform BD simulations using 7 = 0.06 (which is 
a value chosen as a matter of convenience only) and k = 1 as standard values, in a range of 
coupling strengths 20 < F < 1000 covering liquid-to-crystalline states of dusty plasma [i2^ . 

Further, the current density of dust particles is defined as 

1 ^ 

j(r,t) = -^^v,(t)5(r-r,(t)), (2) 

^ i=l 

with the Fourier transform of its cartesian component a given by 

1 ^ 

j„(k,t) = -=^t;,,(t)e^'^-^^W. (3) 

"^^^ i=i 

Assuming that the dust collective modes propagate along the x-direction, i.e., k = (k, 0), 
we let a = X or y, allowing us to define the corresponding longitudinal or transverse current 
density auto-correlation functions as 

Q,tik,t) = {jl{k,0)ji,t{k,t)), (4) 

where the angular brackets denote ensemble averaging over the initial time. The 
longitudinal/transverse spectral densities are then obtained as the Fourier transform of the 
respective current density auto-correlation functions as, 

POO 

Pi,t{k,u)= / dte'^'C^{k,t). (5) 
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The Fourier transforms defined in Eq. ([5]) are evaluated using Fast Fourier transform from 
simulation data. An equivalent and useful definition of the spectral densities, to be used later, 
is given as the real part of the Laplace transform fi271 of the current density auto-correlation 
functions, 

PiAk,^) = 2^{C[Ci,t{k,t)]]s=,^. (6) 
3. Generalized Hydrodynamics 

The basic approach we will take in the development of a GH model for a 2D dusty plasma 
is to treat the dust system as a compressible, viscous neutral fluid of particles interacting via 
the Yukawa potential, analogous to the way how molecules in a Lennard- Jones fluid interact 
ll27l[30l . In that context, we note that the OH approach provides a satisfactory account of the 
long-time (or low-frequency) response of a fluid in the weak coupling regime by means of 
the familiar Naiver-Stokes equation, in which viscous effects are treated as instantaneous (or 
local in time) internal friction force, related to fast thermalization of dust particles due to their 
mutual collisions. However, as the coupling strength increases, the effect of "caging" sets in, 
so that individual dust particles are temporarily trapped in potential wells that migrate through 
the system on a slow time scale, similar to the picture invoked in the QLCA model. Therefore, 
the short-time (or high-frequency) response of the system is dominated by elastic effects due 
to the restoring forces on dust particles in the itinerant potential wells. In the GH approach 
to the Naiver-Stokes equation, a transition from the purely viscous, fluid-like behavior to the 
elastic response of a solid-like system is described by postulating a non-local (in time) friction 
with a memory function characterized by relaxation time r, such that the limit of a viscous 
fluid is recovered at times t ^ r, whereas the solid-like elastic effects are dominant at times 
t < r. 

In addition to the collisions among dust particles, it is necessary to include in the 
GH approach also the effect of their repeated collisions with the neutral molecules in the 
background plasma. These collisions may be treated as a Gaussian white noise, giving rise to 
a picture of Brownian motion, where each dust particle is subject to a local (in time) friction 
force with the neutral drag coefficient 7„. We note that for typical dusty plasmas one expects 
7„r ^ 1. On the other hand, at times t <^ the Brownian motion of dust particles 
due to the neutrals may be considered as ballistic. Therefore, it is justified to ignore the 
effect of neutral drag at the time scale t ~ r -C where viscoelastic effects dominate 
due to interactions among the dust particles [|29l . On the other hand, at long times, such 
that t ~ 7^^ ^ r, both the viscous drag and the neutral drag on dust particles may be 
treated by means of two additive, local in time friction forces in a generalized Navier-Stokes 
equation because the underlying physical mechanisms of energy dissipation are statistically 
independent. 

Therefore, we first ignore the neutral drag and introduce a memory function formalism, 
which incorporates the effects of viscous relaxation into the short-time dynamics of the dust 
collective modes by enforcing high-frequency sum rules upon the power spectral densities 
of such modes [,26 J . The even-order sum rules are defined in terms of the corresponding 
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frequency moments as 

1 

2^ 



duu'^'P^ik^uj) = {uf^ik)), {la) 



whereas all odd moments of frequency vanish since Pi^t{k,uj) are even functions of tu. For 
example, the zeroth-order sum rule is given by [|26ll 
1 

— duPi,tik,u)=vl, (lb) 



271 



— oo 



where vth = ^/kBTd/nid is the dust thermal speed. A connection with microscopic properties 
of the system is accomplished by expressing, e.g., the second moments of frequency for the 
longitudinal and transverse collective modes in Eq. (|7a1 ) with n = 1 in terms of the RDF of 
the dust layer, g{r), respectively as [|26l 

(ufik)) = 3k\t, + £-vl I dhgir) [1 - cos (kx)] (8) 

{oom = k'vtn + d'rg{r) [1 - cos (fcx)] (9) 

Once the initial- value properties of the relevant memory functions are fixed by enforcing 
the sum rules in Eqs. (fTM and (TToI) . one may incorporate the effect of the neutral drag by 
simply adding a local dissipative force into the fluid equations of motion lfT3ll30l . as explained 
in section [ 



3.1. Collective Modes 

The starting point for the extension of OH to GH is given by the linearized continuity equation 
and the Navier-Stokes equation, which may be written for the longitudinal and transverse dust 
current densities as [|26l 

^6pik,t)-tkji{k,t) = 0, (10a) 
ot rUdpXT 

^jtik,t) = -iy,k^jtik,t), (10c) 

where xt is the isothermal compressibility, ui = 2vi + 1J2 is defined as the longitudinal 
viscosity, with ui being related to the shear viscosity and U2 to the bulk viscosity, p is 
the average surface density, 5p is a small perturbation in density, and ji t(k,t) is a small 
perturbation to the longitudinal/transverse current density, assumed to be zero on average. 

Inserting the integral of Eq. (llOal) into Eqs. (|10M and (|10cl) . one can obtain integro- 
differential equations for the current density auto-correlation functions Ci^t{k,t) defined in 
Eq. dH), which must be subject to the initial condition Ci^tik, 0) = vff^. In order to model the 
short-time correlation effects, the resulting equations can be modified to satisfy the frequency 
sum rules in the following way [|25ll26l . In the first step, longitudinal spectral density is forced 
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to satisfy the zeroth frequency sum rule, Eq. (17^ . by replacing the isothermal compressibility 
with the static structure factor, S{k), as 

XT ^4^, (11) 

which generalizes a standard statistical-mechanical relation in the limit /c — i- 0. In the next 
step, generalizations of the longitudinal and shear viscosities are introduced via the as yet 
undefined memory functions, (l)i{k,t) and Kt{k,t), respectively, giving integro-differential 
equations of the form 

ft 



where 



§^C^{k,t) = - ^ dt' K^{k,t-t')C^{k,t'), (12) 

Kiik,t)=^-^^ + k'Uk,t). (13) 
S{k) 

It can be shown that the second frequency sum rule, Eq. (fTal) with n = 1, can be satisfied if 
the initial values of the longitudinal and transverse viscosity memory functions are expressed 
in terms of the corresponding second frequency moments, Eqs. dS]) and respectively, as 

m 

^rA-m-MM ""th .... 



3.2. Ejfects of neutral drag 

For typical conditions in 2D dusty plasmas, the effect of neutral drag on the dispersion 
relations of DAW modes is expected to show up at very low frequencies only, u < 7„, and we 
have indeed found in our BD simulations that the actual value of the drag coefficient 7„ has 
very little effect on the dispersion relations in comparison to the viscoelastic effects. However, 
the spectral widths in Eqs. ([5]) or Q are found to be affected by the neutral drag, as discussed 
in section m 

With the short time dynamics of both the longitudinal and transverse dust collective 
modes fixed by enforcing the sum rules upon the initial-value properties of the relevant 
memory functions, we may now follow suggestions of Refs. lfT3l[30il and introduce the long 
time relaxation effect due to neutral drag by simply adding a local dissipative term in Eq. (fT2)) . 
giving 

§^Ci,t{k, t) = - dt' Ki^tik, t - t')Ci,t{k, t') - 7nG,t(A;, t). (15) 

We proceed with solving Eq. (fT5l) by means of Laplace transform in order to derive 
expressions for the spectral density by invoking the definition Eq. For the longitudinal 
mode we obtain 



' ^ '^[uJ^-Ujl{k)+Ujk''cP'l{k,Uj)Y + {u[k^cP[{k,Uj)+^r]Y' ^ ^ 
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where (t)[{k,uj) and (j)['{k,u) are the real and imaginary parts, respectively, of the Laplace 
transformed longitudinal viscosity memory function, C[(l)i{k,t)]s=iuj, and where we have 
explicitly defined 

Similarly, we obtain for the transverse mode 

with Kl{k, to) and Kl'{k, u) being the real and imaginary parts, respectively, of the Laplace 
transformed transverse viscosity memory function, C[Kt{k, t)]s=i^. 

3.3. Model Memory Functions 

We note that no approximations were used in the development of the GH approach so 
far. Specific forms of Eqs. (fT6l) and (fTSi) may now be obtained by introducing simple 
phenomenological models for the corresponding memory functions. 

3.3.1. Exponential model Choosing an exponential longitudinal viscosity memory function 
of the form (pi{k, t) = (f)i{k, 0)6"*/"^', we obtain fM 



(j)\{k,oo) = Mk,0) ""^ (19a) 

L -\- LU Ti 
2 

^';(k,u) = -Mk,o)-^\^, (1%) 

L + UJ Ti 

to be used in Eq. (fT6l) . We first discuss the limits of very short and very long viscous relaxation 
times Ti. 

In the limit of small relaxation time, ujti -C 1, we recover the so-called delta function 
model for the memory function be letting (f)i(k,0) = ui/ti, where ui is the longitudinal 
viscosity. This situation corresponds to a fluid where viscous drag is described by a local 
friction force, so that Eq. (fT6l) is reduced to ll26l 

Fi[k,uj] = 2v,f, 7T, (20) 

''[u^-uiik)]^ + [uik^ui + ^^)f' 

giving the dispersion relation as w = ujQ{k) with uJo{k) defined in Eq. (flTI) , and having the 
full width at half maximum of 7„ + k'^ui. This result provides a relatively simple account of 
the interplay between the neutral drag and viscous damping, assuming that the quasi-static 
longitudinal viscosity ui can be defined properly ll33ll . 

On the other hand, it is important to study the opposite limit, cuti ^ 1, for the sake of 
comparison with the QLCA model, which is inherently valid in this regime BSl . In particular, 
we find that the viscous damping vanishes in this limit and the elastic effects in the system's 
response prevail, giving the dispersion relation cu = uj^^{k) where 



Dispersion and damping of 2D dust acoustic waves.. 



9 



with the second frequency moment given in terms of the RDF via Eq. ([8]). We note that 
this situation is analogous to that in the QLCA model with two important differences: the 
dispersion relation to = cj^(A;) is different from that found in the QLCA model, and the 
neutral drag still provides a mechanism for damping of the DAW modes. 

For finite values of r;, a dispersion relation in the exponential model for the longitudinal 
viscosity memory function may be obtained from the peak positions in the spectral density 
given by Eq. (fT6l) with Eqs. (|19al ) and (|19Z7| ). This spectral density is fully determined by 
the functions S{k) and {ujf{k)) via Eqs. (fTTI) and (I14al) . respectively, and by the longitudinal 
viscosity relaxation time ti, which may be a function of k as well. We note that the former 
two functions can be calculated, at least in principle, from the RDF of the dust layer by using 
the usual definition for S{k), and Eq. dS]) for {uf{k)). However, while g{r) can be obtained 
directly from BD simulations, or may be available from first principles, there is no simple way 
to determine ti from first principles. One possibility to proceed is to ignore its k dependence 
and treat r; as a free parameter, possibly dependent on F, that can be determined from, e.g., 
fitting the peak positions of Eq. (fT6l) to the peak positions of the spectral densities obtained 
from experiments or computer simulations. It is expected that such a procedure would yield 
a dispersion relation lying somewhere between u = uio{k) and u = uo^^{k), defined via Eqs. 
(fTTI) and (I2TI) . respectively. 

Implementation of the exponential model for transverse viscosity memory function 
Kt{k,t) is a straightforward repetition of the procedure outlined above, and it produces 
a result for the spectral density that is entirely equivalent to that derived by Murillo for 
transverse modes in 3D strongly coupled Yukawa liquids (T^. In brief, one assumes 
Kt{k, t) = Kt{k, 0)6"*/"^' with the initial value given in Eq. (I14ZpI) . where the second moment 
of frequency for the transverse mode may also be obtained from the RDF by using Eq. dH). 
One further obtains the real and imaginary parts of the Laplace transform C[Kt{k, t)]s=j^ by 
using expressions analogous to Eqs. (|19al) and (|19Z7|) . which, upon substitution into Eq. (fTSl) 
yield a spectral density for the transverse mode in the exponential model. 

It should be reiterated that, in the case of vanishing relaxation time, = 0, the transverse 
mode is purely diffusive and the resulting spectral density in Eq. (fTSi) exhibits no dispersion 
ll26l . In the opposite limit of — )• 00, viscous damping vanishes, while the dispersion 
relation, given by w = with 



JJk) = vW. (22) 

Vth 

can be shown to exhibit a quasi-acoustic behavior in the long wavelength limit. On the other 
hand, when using the peak positions of the spectral density in Eq. (fTSi) for the exponential 
model with finite r^, one encounters a cutoff wavenumber, kc, in the dispersion relation for 
the transverse mode, which can be estimated from the condition {ujf{kc)) = {vth/TtY in the 
limit of vanishing neutral drag [fT4ll22| . However, like in the case of the longitudinal mode, 
since there is no simple way to determine from first principles, one can again treat it as a 
free parameter that may be determined from an appropriate fitting, e.g., by using the cutoff 
values kc observed in the experiments on shear waves in 2D Yukawa liquids [[HI. 
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3.3.2. Gaussian model Instead of using r; and as free parameters, one may attempt 
enforcing higher-order frequency sum rules to see if a refinement of the above exponential 
model can be achieved. However, since the exponential model for the viscosity memory 
functions does not support moments higher than the second order when used in Eqs. (fT6l) 
or (fTSi) . one needs a model with different time dependence. Besides having to satisfy the 
second-frequency sum rule via Eq. (|14al) . it is shown in the Appendix that such model must 
additionally satisfy both the third-order sum rule giving \^^^4>i{k, = 0, and the fourth- 

order sum rule giving 



"th 



J t=0 ^'^'^th 



(23) 



We see now that a Gaussian model for the longitudinal viscosity memory function of 
the form t) = (pi{k, 0) exp [—t'^/crfik)] will satisfy all moments up to and including the 
fourth if the Gaussian relaxation time cri(k) is given by 

with the initial value of the memory function, </);(A;,0), given by Eq. (I14al) . Finally, by 
using the Laplace transform of the Gaussian longitudinal viscosity memory function with 
s = icu in Eq. (fT6l) . we obtain a spectral density which is fully determined by three functions: 
S{k), {ujf{k)), and {ujf{k)) via Eqs. (flTI) and (|14al) without the need for free parameters. 
Unfortunately, analytical expressions that can be used for calculation of the fourth frequency 
moment are rather cumbersome and require a three-particle distribution function, which is 
difficult to obtain from simulations [261 . Therefore, since calculation of {ujf{k)) from first 
principles is impractical, we may use the Gaussian model by computing the fourth moment 
numerically from Eq. (iTal) with n = 2 where Pi{k,uj) is obtained from simulation. For 
consistency, it is then desirable to also compute the second moment from Eq. (iTal) with n = 1 
using the same Pi{k,uj) based on simulation data. In this way, we may use the Gaussian 
model as a simulation-based, parameter- free test of the quality of the approximation achieved 
in using the exponential model with a suitable choice of finite relaxation time r^. 

We finally note that a completely analogous development of the Gaussian model can be 
applied to the transverse mode. 



4. Results and Discussion 



We tested several values for the reduced neutral drag coefficient 7 in our BD simulation and 
found no noticeable dependence of the resulting spectra on 7. Moreover, all the quantities 
used in the GH modeling of the collective modes (e.g., the RDF and the relaxation times in 
the exponential model of the memory functions for the longitudinal and transverse modes) 
were also found to be robustly independent of the (small) values of 7 used in simulation. 
Therefore, all results will be shown for a standard value of 7 = 0.06, along with the standard 
screening parameter k = 1. 
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Figure 1. Radial distribution function g{r) as a function of the reduced distance r/a 
for K = 1.0, 7 = 0.06, and T = 20, 60, 100, 200, 600, and 1000, corresponding to the 
increasing hight of the principal peak, respectively. 



4.1. longitudinal Wave Mode 

We use coupling strength with values F = 20, 60, 100, 200, 600, and 1000 to investigate the 
longitudinal mode in a broad range of dusty plasma conditions, going from liquid to crystalline 
states. 

As described above, both the static structure factor S{k) and the second frequency 
moments can be calculated from the equilibrium RDF g{r), which is shown in Fig. [Hfor 
several values of F. However, there are significant difficulties related to using S{k) for 
modeling the dispersion relation of the longitudinal mode at long wavelengths. Namely, 
when S{k;) is calculated from the RDF, its small-A; values exhibit a very strong dependence 
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r=20 
r=60 
r=ioo 




Figure 2. Static structure factor S{k) as a function of the reduced wave-number ka for 
K = 1.0, 7 = 0.06, and T = 20, 60, 100, 200, 600, and 1000, corresponding to the 
increasing hight of the principal peak, respectively. 



on the fluctuations appearing in the simulation due to the finite number of dust particles 
[|34l . This is expected to give rise to rather noisy dispersion curves, especially in the 
situations characterized by short relaxation times when the dispersion is dominated by the 
frequency ijjQ{k) defined in Eq. (fTTI) . To remedy the situation to some extent, we resort 
to calculating S{k) directly from the simulation data by using the long time average of the 
Fourier transformed particle density in equilibrium, as follows 

S{k) = (p*(k)p(k)), (25) 
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Figure 3. Longitudinal wave dispersion curves against simulation spectra for k = 1.0, 
7 = 0.06, and T = 20, 60, 100, 200, 600, and 1000. The upper (dark green), middle 
(red) and the lower (light green) solid curves represent the exponential model with 
an infinitely long, finite, and vanishing relaxation times, respectively. The thick solid 
(black) curve represents the Gaussian model for wave-numbers ka < 6. 



where 



P(k) 



AT 



,ikr,: 



(26) 



In Fig.[2l we show the thus obtained results for S{k) for several T values, and note that the 
peak structures in that figure were found to be identical to those arising in S{k) computed 
from the corresponding RDFs in Fig.fT] 

In Fig. |3] we show the simulation data for the spectral density of the longitudinal mode, 
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along with the dispersion curves from the exponential model, represented by three solid curves 
corresponding to, in descending order, infinitely long, finite, and vanishing relaxation times 
T/. Also shown in Fig. [3]by a thick black solid line is the dispersion curve obtained from the 
Gaussian model with the second and fourth frequency moments evaluated from Eq. (ITal) with 
n = 1 and 2, respectively, where Pi{k,uj) was taken directly from the simulation data for 
spectral density. While the three curves for the exponential model extend over the full range 
of wave-numbers shown in Fig. [31 the results for the Gaussian model are shown over a reduced 
range of wave-numbers covering the first Brillouin zone only, ka < 6, because computations 
of the frequency moments via Eq. (iTal) were hampered by a significant noise in the simulation 
spectra at larger wave-numbers and lower coupling strengths. 

One notices in Fig. [3] that, for the exponential model with vanishing relaxation time 
(or, equivalently, the delta- function model), the resulting dispersion relation, cu = ujo{k) with 
ujo{k) given in Eq. (flTl) . displays some noise coming from the simulation data via Eq. (l25l) . We 
found this noise to be significantly weaker than the noise arising when S{k) is calculated from 
the RDF. On the other hand, as the longitudinal relaxation time ti increases in the exponential 
model, there is a relative increase in the contribution of the second frequency moment {ujf{k)), 
which, when evaluated from the RDF via Eq. ([8]), exhibits a much smoother k dependence than 
uio{k) with S{k) evaluated from Eq. (l25l) . 

As mentioned before, the dispersion relation for finite values of the relaxation time ti 
in the exponential model for the longitudinal viscosity memory function is obtained from 
the peak positions in the spectral density given by Eq. (fT6l) with Eqs. (|19al ), (|1%I) . (flTI) . 
and (I14al) . We have chosen the values of ti which provide the best fit to the peaks in the 
simulation spectral densities, and we have found ti to be a relatively weak function of T that 
may be reasonably well approximated by 



We emphasize that the quality of our choice of the values for ti in the exponential model is 
confirmed through a close agreement with the dispersion curves from the Gaussian model for 
wave-numbers in the first Brillouin zone, as displayed in Fig.[3l 

It is worth mentioning that the extreme cases of the exponential model, corresponding 
to the limits — )• and ti — )• oo, yield two parameter-free dispersion relations for the 
longitudinal mode, to = ujo{k) and tu = uj^^{k), respectively. Remarkably, it is noticed in 
Fig. [3] that these two dispersion relations provide good account of the simulation data by 
straddling the thermal noise, seen in the recorded spectra at the low-to-medium V values. 
Moreover, at the medium-to-high V values, one notices that u = uJo{k) closely follows the 
dispersion relation from the Gaussian model for the wave-numbers in the first Brillouin zone 
and, in particular, reproduces the near-vanishing of frequency at around /ca ~ 4 for F = 600 
and 1000. This is somewhat surprising, given that the delta-function model is stretched well 
into the condensed state at those two coupling strengths. 

We further compare in Fig.|4]the dispersion relations u = uJo{k) and u = ujlo{k) (shown 




for r < 60 



for 200 < r < 1000. 



for 60 < r < 200 



(27) 
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Figure 4. Longitudinal wave dispersion curves against simulation spectra for k = 1.0, 
7 = 0.06, and T = 20, 60, 100, 200, 600, and 1000. The upper (black) and 
the lower (red) dashed curve represent the results from the EQLCA and QLCA 
models, respectively. The upper (dark green) and the lower (light green) solid curves 
represent the exponential model with an infinitely long and vanishing relaxation times, 
respectively. 
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by the upper and lower solid curves, respectively) with the results from the QLCA model 
ll24l . shown by the lower dashed curve. It is immediately obvious that the QLCA dispersion 
does not reproduce the direct thermal effect, which is responsible for a quasi-linear increase 
in the peak frequencies of the simulation spectra at higher wave-numbers and lower coupling 
strengths. In a previous work, we have shown that this deficiency can be easily rectified by 
a simple extension of the QLCA model, labeled as the EQLCA model in Ref.[22i], which 
is shown by the upper dashed curve in Fig. IH One notices a surprisingly good agreement 
between the EQLCA dispersion relation and the curve uj = colo{k) from the exponential 
model for — )• oo. While an agreement with the QLCA dispersion relation was expected 
at the higher T values and lower k values owing to the fact that the QLCA model inherently 
assumes cuti ^ 1 [1351 and the direct thermal effect is relatively weak at low k and high T 
values [|22||. it is remarkable how the parameter- free dispersion relation co = uj^^{k) provides 
justification for the empirically derived EQLCA model over the broad ranges of wavenumbers 
and coupling strengths. 

We next compare the performances of the exponential model with finite relaxation time 
Ti from Eq. (l27l) and the Gaussian model with the frequency moments evaluated from the 
simulation spectra via Eq. (fTal) by using these two models in Eq. (fT6l) to predict the frequency 
dependent spectral profiles at several fixed wave-numbers. Given that the half-width at half- 
maximum (HWHM) of these profiles may be directly related to the wave-number dependent 
damping rate of the longitudinal DAW mode, in this way we demonstrate the advantage of 
using the memory function formalism within the GH model in tackling the difficult issue of 
damping of the collective excitation modes in dusty plasmas. 

Specifically, we show in Figs. 5-10 the frequency dependencies obtained from the 
spectral density in Eq. (fT6l) with the exponential and the Gaussian models for four wave- 
numbers, ka = 0.61, 1.17, 2.29, and 3.36, and compare them with the corresponding profiles 
of the simulation spectral density for T = 20, 60, 100, 200, 600, and 1000, respectively. 
One notices in Figs. 5-10 a fair agreement between the simulation data and both GH models, 
which is especially good for lower F values where the hydrodynamic regime is expected to be 
more pronounced, while for very high coupling strengths, say, F > 200, the agreement begins 
to deteriorate. Comparison of both GH models with the simulation data may be considered 
fairly reasonable, even at high coupling strengths where the hydrodynamic model has been 
stretched well into the crystalline state. 

Finally, in Fig.[TTl we show a comparison between the wave-number dependent HWHM 
values, obtained from the simulation spectral density profiles and from Eq. (fT6l) where we 
used both the exponential model with finite ti values given in Eq. (l27l) . and the Gaussian 
model with the frequency moments evaluated from the simulation spectra via Eq. (l7al) . One 
notices that all three sets of data are noisy due to the noise in the simulation spectra (one 
recalls that, in the case of the exponential model, the noise stems from ujo{k) with S{k) 
evaluated from Eq. (|25l) ). As expected, the HWHM values from the Gaussian model are 
systematically somewhat larger than those from the exponential model, but they are both seen 
to be in reasonably good agreement with the simulation data, even though the agreement 
becomes blurred by the increased noise at the highest F values, say F > 200. One also notices 
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Figure 5. Longitudinal spectral density profile curves versus reduced frequency oj / ujpd 
for K = 1.0, 7 = 0.06, r = 20, and ka = 0.61, 1.17, 2.29, and 3.36. Simulation data 
are shown by the noisy (blue) solid curve, exponential model with finite relaxation 
time is represented by the dashed (red) curve, and the Gaussian model is represented 
by the smooth solid (green) curve. 



in Fig.[TT]that all damping rates roughly follow a quadratic dependence on k, with increasing 
opening of the parabola as T gama increases, and with the limiting value as /c — )■ being close 
to the damping rate due to the neutral drag, given by 7/2 = 0.03. 



4.2. Transverse Wave Mode 

Transverse, or shear waves are expected to only exist for higher coupling strengths T [fT4l [8l. 
In the following we compare the simulation results with the GH approach for a typical case 
of r = 100. 
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Figure 6. The same as Fig. |5l but for T = 60. 



Using Eq. ^ to evaluate the second frequency moment from the RDF, we have obtained 
a dispersion relation from the peak positions of the spectral density in Eq. (fTSi) for the 
exponential model with = 5.0 (note that typical values of Tt for the transverse mode 
are generally found to be larger than those for the longitudinal mode, in agreement with the 
remarks of Boon and Yip ^26^). The results are shown in Fig. [121 along with the dispersion 
relation in the limit of infinite relaxation time, — t- oo, given by w = w^(fc) with Eq. (|22|) . 
and are compared with the simulation spectrum for V = 100. Both theoretical dispersion 
curves show good agreement with the simulation data for ka < A and, while the simulation 
spectrum is not conclusive about the cutoff wave-number, the model with finite Tj clearly 
exhibits a cutoff at fcc 0.5/a [Hid. 

Further, several profile plots of the spectral density are shown in Fig. [13] for F = 100, 



Dispersion and damping of 2D dust acoustic waves.. 



19 




Figure 7. The same as Fig.|5l but for T = 100. 



where the results from Eq. (fTSi) for both the exponential model with = 5.0 and the Gaussian 
model with the frequency moments evaluated from the simulation data via Eq. (TTal) are 
compared with the simulation profiles for fca = 0.61, 1.17, 2.29, and 3.36. One notices that the 
agreement of the exponential model with the simulation spectral profiles is quite satisfactory 
for the lower two wave-numbers, but the simulation profiles at the two higher wave-numbers 
are seen to be much broader than the corresponding peaks from the exponential model. On 
the other hand, the Gaussian model does not reproduce the dispersion of the transverse mode 
at all, since the corresponding profile curves peak at zero frequency for all wave-numbers, as 
displayed in Fig. [131 Nevertheless, the broadening of the simulation profiles at the two higher 
wave-numbers in Fig. [13] seems to be better echoed by the Gaussian than by the exponential 
spectral profile curves. 



Dispersion and damping of 2D dust acoustic waves.. 



20 




Figure 8. The same as Fig. |5l but for T = 200. 



5. Concluding Remarks 

We have carried out Brownian Dynamics simulation of a 2D layer of strongly-coupled dusty 
plasma and extracted the equilibrium radial distribution function, static structure factor, and 
the spectral densities for both the longitudinal and transverse dust acoustic wave modes. 

We have then shown that the memory function formalism of the theory of generalized 
hydrodynamics provides good semi-analytic results to model these wave modes. In particular, 
following the approach of Boon and Yip [|26l . we have developed an exponential model for the 
viscosity memory functions that satisfies the second frequency sum rule with the relaxation 
time being used as a fitting parameter. With such exponential model we have obtained 
good fits to the simulation data for the longitudinal dispersion curves over a wide range of 
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Figure 9. The same as Fig. |5l but for T = 600. 



coupling strengths, 20 < F < 1000, as well as reasonable match with the spectral density 
profiles for several values of the wave-number k. Next, we have extended the theory in 
a straightforward manner to satisfy the third and fourth frequency sum rules, providing a 
parameter-free Gaussian model for the viscosity memory function. The results using this 
model provided a successful test for our choice of the relaxation time in the exponential model 
in modeling both the dispersion relations and the spectral density profiles. 

Moreover, we pointed out that the limits of an infinitely long and infinitesimally short 
relaxation times in the exponential model represent two parameter-free models that are fully 
determined by the equilibrium radial distribution function and the static structure factor of the 
system. It was then shown that these two limits provide good upper and lower bounds for 
the thermal noise observed in the simulation spectra for the longitudinal mode. Comparisons 



Dispersion and damping of 2D dust acoustic waves.. 



22 




Figure 10. The same as Fig.|5l but for T = 1000. 



of the dispersion relations obtained from these two limits of the exponential model with the 
dispersion relations obtained from both the Quasi-localized charge approximation (QLCA) 
and the extended QLCA (EQLCA) showed that the limit of an infinitely long relaxation time 
reproduces remarkably well all the features of the EQLCA model, including the direct thermal 
effect. On the other hand, the limit of an infinitesimally short relaxation time was found to 
reproduce the near- vanishing of the dispersion relation, which is seen in the simulation spectra 
at certain short wave-lengths for large coupling strengths, but is not reproduced by the QLCA 
group of results. 

In addition, we have also tested the memory function formalism within the theory of 
generalized hydrodynamics against the simulation results obtained for transverse wave modes 
at the coupling strength of F = 100, and found that the exponential model with relaxation time 
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Figure 11. Half-width at half maximum of the spectral profile curves versus the reduced wave- 
number ka for K. = 1.0, 7 = 0.06, and T = 20, 60, 100, 200, 600, and 1000, with the (blue) 
dots representing the simulation data, the (red) noisy solid curve representing the exponential 
model with finite relaxation time, and the squares representing the Gaussian model. 



used as a fitting parameter can provide a reasonable estimate for the wave dispersion relation, 
including the appearance of a cutoff wave-number. Both the peak positions and the widths 
in the spectral density profiles from the simulation were well reproduced by the exponential 
model at long wavelengths, whereas the broadening of these spectra at short wavelengths was 
better echoed by the Gaussian model. 

The most important finding of the present analysis is that the wave-number dependent 
damping rates, which were extracted from the simulation spectra, were found in good 
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Figure 12. Transverse wave dispersion curves against the simulation spectra for 
K = 1.0, 7 = 0.06, and F = 100, with the upper solid (green) curve and the lower 
solid (blue) curve representing the exponential model with an infinitely long and finite 
relaxation times, respectively. 

agreement with the results from both the exponential and Gaussian models, testifying to 
the strengths and potentials of using the memory-function formalism within the GH theory 
in describing the damping of the longitudinal dust acoustic wave in strongly-coupled dusty 
plasmas. 

Finally, we have included the effect of neutral drag due to the collisions of dust particle 
with neutral molecules in the background plasma via a (local in time) drag coefficient. While 
we did not find that the neutral drag affects our modeling of dispersion relations for both the 
longitudinal and transverse waves in any significant way, the overall damping rates for the 
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Figure 13. Transverse spectral density profile curves versus reduced frequency a; / ujpd 
for K = 1.0, 7 = 0.06, r = 100, and ka = 0.61, 1.17, 2.29, and 3.36. Simulation data 
are shown by the noisy (blue) solid curve, exponential model with finite relaxation 
time is represented by the dashed (red) curve, and the Gaussian model is represented 
by the smooth solid (black) curve. 



longitudinal wave showed possibly important roles of the neutral drag at long wavelengths 
and high coupling strengths where viscous damping is expected to be reduced. This aspect of 
the GH theory shall be studied in future work. 

In conclusion, the memory function formalism of the generalized hydrodynamics 
provides a very promising route to analytical modeling of spectral density of collective modes 
in strongly coupled 2D dusty plasmas and, in particular it demonstrates a good account of 
viscoelastic effects which are not well described by other analytical methods. 
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Appendix A. Frequency sum rules 



Initial- value properties of the current density auto-correlation functions are related to the even- 
order frequency moments via ll26l 

Q2n 



dt 



2n 



CUk,t) 



(A.l) 



i=0 



whereas all odd-order moments must vanish. We extend here the development outlined by 
Boon and Yip [12611 by introducing the third and fourth frequency sum rules into the GH theory 
of spectral density. Since these sum rules are used to fix the short time properties of the model 
memory functions, we may again initially neglect the neutral drag. 

To illustrate this procedure for the longitudinal mode, we differentiate Eq. (fT2l) twice to 
obtain the third derivative as 



^3 



d 



e j\t'Ci{kX)^Hk,t-t') 



Ci{k,t')-Mk,t-t') 



J t'=t 



e^i(k,0)^Qik,t) 



which on evaluating at t = 0, yields 



Q{k,t) 



t=0 



d 



Ci{k,t')-Mk,t-t') 



t'=t=o 



Since the third moment must vanish, we have the condition 



k'Ciik,0) 



^Mk,t-t') 



0, 



(A.2) 



(A.3) 



(A.4) 



t'=t=o 



showing that the first derivative of the longitudinal viscosity memory function must vanish at 
initial time. This is to be expected because a memory function should be viewed as being just 
another time auto-correlation function with vanishing derivatives of odd order at the initial 
time ll26l . Next, evaluating the fourth moment we find 



s{k) dt^ 



r (9^ 

Clik, t)-k^ J dt' Qik, t')^Mk, t - t') 



k' 



rj2 

Q(k,t')—Mk,t-t') 
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dt 
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dt^ 

Again setting t = we obtain. 



(A.5) 



dt^ 
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(A.6) 



t=0 
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Substituting 0;(A;, 0) from Eq. (|14a|) and invoking Eq. (lA.lD with n = 1 and 2, we obtain Eq. 
(l23l) in the main text. 
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